#Replication Material for Generalizability of IR experiments beyond the U.S.
#Produce figure for prospective and retrospective power analyses

library("tidyverse")
library("broom")
library("estimatr")
library("ggplot2")
library("metafor")
library("ggthemes")
library("ggh4x")
library("haven")
library("ggimage")
library("devtools")
install_github("jimjam-slam/ggflags", dependencies = TRUE)
library("ggflags")
library("texreg")
library("ltm")
library("interflex")
library("stargazer")
library("tidytext")
library("wordcloud")
library("gridExtra")
library("quanteda")
install_github("naoki-egami/exr", dependencies = TRUE)
install_github("naoki-egami/evalid", dependencies = TRUE)
library("evalid")
library("exr")
library("hettx")
library("DeclareDesign")
library(tidyverse)
library(ggthemes)
set.seed(999)

# Simulate Data -----

design <- 
  declare_model(
    country = add_level(N = 7),# Declare 7 countries. 
    subjects = add_level(N = 2800, # Declare 2800 participants per counry
                         U = rnorm(N), # random error
                         potential_outcomes(Y ~ -0.12*Z + U)) 
  ) +
  declare_assignment(Z = block_ra(country)) + # Block randomize treatment by country
  declare_measurement(Y = reveal_outcomes(Y ~ Z)) +
  declare_estimator(handler = function(data) { # Declare estimator, country specific regression.
    data |>
      group_by(country) |>
      summarize(tidy(lm_robust(Y ~ Z)), .groups = 'drop') |>
      filter(term == "Z")
  })

simulations <- simulate_design(design, sims = 500)


# Prospective analysis -----

# fraction of sims reject null with adjusted p
simulation_pcorrections <-
  simulations %>% 
  group_by(sim_ID) %>% 
  mutate(.,
         adjust_p_bh = p.adjust(p.value, method = "BH"),# Implement BH correction per simulation
         adjust_p_bnf = p.adjust(p.value, method = "bonferroni"),# Implement Bonferroni correction per simulation
         sig_exp = ifelse(p.value <0.05,1,0),# Create naive indicator of hypothesis rejection
         sig_exp_adj_bh = ifelse(adjust_p_bh <0.05,1,0),# Create indicator of hypothesis rejection, with BH
         sig_exp_adj_bnf = ifelse(adjust_p_bnf <0.05,1,0),# Create indicator of hypothesis rejection, with Bnf
  ) %>% 
  ungroup()


mean(simulation_pcorrections$sig_exp) # Power unadjusted p values 0.88
mean(simulation_pcorrections$sig_exp_adj_bh)# Power BH p values 0.89
mean(simulation_pcorrections$sig_exp_adj_bnf)# Power Bnfrni p values 0.69


# Plot distribution of simulation p values
ggplot(simulation_pcorrections, aes(x=adjust_p_bh)) + 
  geom_histogram()+
  geom_vline(aes(xintercept=0.05),
             color="red", linetype="dashed", size=1)+
  labs(x="BH Adjusted P Value",
       y="Number of Simulations")+
  annotate("text", x = 0.75, y = 2000, label = "500 Simulations \nHypothesized ATE: -0.12 (IL ATE) \n 2800 Subjects per country (7 countries) \nproportion p (adjusted BH) < 0.05 = 0.88 \nproportion p (non-adjusted) < 0.05 = 0.89 \nproportion p (adjusted Bonferroni) < 0.05 = 0.69") +
  theme_bw()
ggsave("../figures/appendix/power.pdf",width = 8, height = 4)


# Implement Sign Generalization Test for simulation (retrospective analysis) -----

simulations<- simulations %>% 
  mutate(., z = estimate/std.error,
         dem_sign_gen = pnorm(z))

  
merged <- data.frame()
for(x in 1:500){
a <-simulations %>% 
    filter(sim_ID==x)

b<-pct(a$dem_sign_gen)

merged<-bind_rows(merged,b)
}


sim_frac <-
  merged %>%
  mutate(.,
         sig_exp = ifelse(p_value <0.05,1,0),# Create naive indicator of hypothesis rejection
  ) %>% 
  ungroup()


mean(sim_frac$sig_exp) # Power p values 0.9

# Plot distribution of simulation p values
ggplot(merged, aes(x=p_value)) + 
  geom_histogram()+
  geom_vline(aes(xintercept=0.05),
             color="red", linetype="dashed", size=1)+
  labs(x="Sign Generalization P Value",
       y="Number of Simulations")+
  annotate("text", x = 0.4, y = 2000, label = "3,500 simulated studies \n(7 countries*500 simulations) \nHypothesized ATE: -0.12 \n2,800 subjects in each country simulation \nproportion p (partial conjunction) < 0.05 = 0.9") +
  theme_bw()
ggsave("../figures/appendix/power_simulation.pdf", width = 8, height = 4)

